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Abstract 

We present a new method to calculate directly the one-particle self-energy of an impurity 
Anderson model with Wilson's numerical Renormalization Group method by writing this 
quantity as the ratio of two correlation functions. This way of calculating T,(z) turns 
out to be considerably more reliable and accurate than via the impurity Green's function 
alone. We show results for the self-energy for the case of a constant coupling between 
impurity and conduction band (9raA(w + i0 + ) = const) and the effective A(z) arising in 
the Dynamical Mean Field Theory of the Hubbard model. Implications to the problem 
of the metal- insulator transition in the Hubbard model are also discussed. 

1 Introduction 

The single impurity Anderson model [jT| is one of the most fundamental and probably the 
best understood model for strong electronic correlations. Invented to describe the properties 
of magnetic impurities in non-magnetic metallic hosts 35 years ago, a variety of standard 
techniques have been applied to it and new methods have been developed to study its static and 
dynamic properties in basically the whole parameter space (for a review see e.g. @). Although a 
very clear picture of the physics of the single impurity Anderson model has emerged from these 
calculations, a reliable method for calculating dynamic properties at very low temperatures and 
intermediate or large values of the Coulomb interaction was for a long time missing. 

For example, Bethe ansatz calculations ||, which are essentially exact, can only access 
static properties and the Quantum Monte Carlo method JO, which can be viewed as another 
numerically exact technique, cannot reach very low temperatures and/or large values of the 
Coulomb parameter, although it does not suffer from a minus sign problem here. In addition, 
the analytic continuation of the imaginary time data to real frequencies is a numerically highly 
ill-conditioned problem. 

Among the approximate treatments the resolvent perturbation theory together with the 
so-called Non-Crossing Approximation || turned out to be a simple and powerful technique 
for high and intermediate temperatures of the order of the Kondo scale but completely fails 
to reproduce the local Fermi-liquid properties as T — > 0. Last but not least, straightforward 



1 



second-order perturbation theory in U @ has been shown to work surprisingly well down to 
T = but it is restricted to the symmetric case and not too large values of U. 

The Numerical Renormalization Group method (NRG), invented by Wilson for the Kondo 
problem |7J and later applied by Krishnamurthy et al. to the impurity Anderson model ||, 
is usually acknowledged primarily in the context of universality and low-energy fixed point 
behaviour of the Kondo or Anderson model. One of its most appealing features is that it can 
deal equally well with small, intermediate or large values of U and is not restricted to half- 
filling. During the last 15 years considerable progress has been made to extract dynamical 
properties with this method, too, and it has been shown to give very accurate results also for 
e.g. dynamical one- and two-particle and also transport properties |9|, [T(J. The NRG works 



best at T = 0, and various dynamic correlation functions can be calculated with an accuracy 
of a few percent. Although less well defined for finite temperatures, its extension to T>0 also 
shows very good agreement with exact results [ID|. It is quite remarkable that no sum-rules 



(Friedel sum rule, total spectral weight) must be used as input for these calculations. On the 
contrary, they can serve as an independent check on the quality of the results. 

More recent interest in reliable methods to solve the impurity Anderson model and calculate 
its dynamical properties has been motivated by the discovery that lattice models in the limit 



of infinite dimensions acquire a purely local one particle self-energy |]TTj. This simplification 
eventually leads to a mapping of the lattice problem on an effective impurity Anderson model 
coupled to a medium to be determined self-consistently ||l2 |. Note that in the general case 



the solution of this self-consistency requires the knowledge of the one particle self-energy. In 
view of the wide range of problems to which this so-called Dynamical Mean Field Theory 
(DMFT, see e.g. ||13|| ) can be applied, it seems surprising why there have been hardly any 
contributions using the NRG. The only NRG-calculation known to us is the work of Sakai et al. 
14 1 where the symmetric Hubbard model in the metallic regime was studied. In their paper, 



these authors point out some difficulties in the progress of iterating the NRG results with the 
DMFT equations, which are largely related to the necessary broadening of the NRG spectra 
(see further below). 

In this contribution we present a new method to calculate dynamic properties for the im- 
purity Anderson model, namely by directly constructing the interaction contribution to the 
self-energy as the ratio of two correlation functions, X^(z) = UF a (z)/G a (z), with F a (z) = 
((faf^fa, fa))z and G a (z) = ((f a , fl)) z (see Section II). Details of how to calculate the F(z) are 
given in the appendix. In section III we discuss results for 

• the standard case, where the coupling between impurity states and metallic host, QmA(uj+ 
i0 + ), is constant, 

• and the Hubbard model in d = oo, where A(z) has to be determined self-consistently. 

The Hubbard model is studied in the paramagnetic regime, at half-filling and T = 0. We discuss 
the properties of self-energy and local density of states both in the metallic and insulating 
regimes and some preliminary results for the metal-insulator transition. 
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2 Calculation of the self-energy 



Model and basic concepts 

The impurity Anderson model is written in the form 

(7 

+ J2 £ k C L C k«+Y. V k{fi C k« + cLf«)- (i) 

ka ka 

In the model ([!]), denote standard annihilation (creation) operators for band states with 
spin a and energy e*., those for impurity states with spin a and energy £f. The Coulomb 
interaction for two electrons at the impurity site is given by U and both subsystems are coupled 
via a hybridization 14, which we allow to be /c-dependent here. 

Our final goal is to calculate the one-particle Green's function G a (z) = ((fa,f a ))z, which 
formally can be written as 

W = • < 2 > 

While this formal introduction of the one particle self-energy £(V) is straight forward, the actual 
calculation of G(z) or alternatively £(z) usually is an extremely complicated problem. In order 
to express the self-energy £(z) by standard impurity correlation functions, we make use of the 
equation of motion 

z((A,B)) t + ((CA,B)) z = ([A,B] v ), (3) 

with £■ = [H, •]_ and r] = +, if both A and B are fermionic operators, 77 = — otherwise. The 
correlation functions are defined as ((A, B)) z = i J °° e tzt ([A(t), B] v ). For A = f a and B = /J we 
obtain the equation of motion for the f-Green's function as 

(z - e { )G a (z) - l/fl/jj/,, fl ))z _ J- y fc((cfcCT; /t^ = x . (4) 

The correlation function ((c). a , f£)) z is related to G a (z) via eq. (^|) with A = c^ a and B = /t 
through 

(^-^((Wj>>,- W CT W = 0. (5) 

The £7-term does not enter this equation as the Coulomb interaction only acts on the impurity 
states. Together with ([|) eq. @) has the form 

(z - s { )G a (z) - UF a (z) - A(z)G t7 (z) = 1, (6) 

where we have defined F a (z) = ((fafafa, f a ))z and A(z) = Hk^k^~- The total self-energy 
Eo^z) for the single impurity Anderson model is thus given by 

V a {z)=A(z) + I%(z) , (7) 
where the nontrivial part due to the Coulomb correlations (z) is obtained from 

=?W = V§£ (8) 
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For simplicity and since we are only interested in the paramagnetic situation for the time being 
the spin index will be dropped in the following. 

Alternatively, the interaction part of the self-energy can of course also be calculated directly 
from eq. using 

YF{z) = Go(z)- 1 - G{z)-\ with G (z) = 1 —u^- (9) 

On a first glance, there seems to be no apparent reason to prefer the more complicated equation 
(|8|) over equation (|9|). In order to clarify the advantage of using eq. (|j) instead of eq. (Jffy for 
the calculation of TP '(z) with the NRG we want to give a brief description of how the spectral 
densities for G(z) and F(z) are calculated with the NRG. 



Technical details 

Within the NRG, the impurity Anderson model eq. ([!]) is mapped onto a semi-infinite 
chain (see |7|, ||) which is diagonalized iteratively starting from the uncoupled impurity. At 
each iteration, the number of states increases by a factor of 4 and after a certain number of 
iterations, the basis kept for the next iteration has to be truncated. The important point of the 
method is that the coupling between consecutive elements of the chain decreases exponentially 
for increasing distance from the origin, so that with increasing chain length at each iteration 
basically only the lowest lying states will be renormalized and such a truncation is meaningful. 
The spectral functions at each iteration are calculated from the corresponding matrix elements, 
which are in turn related to those of the previous iteration. This procedure is well established 
for the one-particle density of states A(u) = —- ( 3mG(uj+iO + ) || [lOj and can straight forwardly 
be extended to B{uj) = —^mF(u + i0 + ). For details we refer the reader to the appendix. 
Due to the truncation of states, the spectral function for the whole frequency range has to be 
built up from the data of all the iterations. 

The resulting spectral function is a set of ^-functions at frequencies uo n with weights g n 
which are broadened on a logarithmic scale as 



g n 5(uj - u n ) 



9n 



n 



exp 



(In to — In uj v 



(10) 



This form of broadening was also used in || and Jl0| and is especially adapted to the exponential 
variation in energies peculiar to the NRG. The width b n is chosen as b independent of n and 
we use values 0.3 < b < 0.6. 

It is well known that with this scheme the NRG gives already quite accurate results for G(z) 
H p!T5[1 . However, one might anticipate some problems with the calculation of TP (z) using eq. 
(|^). The function G (^) _1 is, of course, known exactly since A(z) is a given quantity. Building 
the difference between an exactly known and a numerically determined function is usually very 
susceptible to numerical errors, especially in regions where the result becomes small. Since this 
is expected to happen close to the Fermi level, i.e. in the physically most relevant region, one 
is likely to run into problems there. 

One naive attempt to reduce these kind of inconsistencies and numerical errors when build- 
ing the difference in eq. @ is to treat Gq(z)~ 1 and G{z)~ 1 on the same level, that is to calculate 
Gq{z)~ 1 via the NRG as well by setting U = 0. However, since according to the theory of error 
propagation in sums or differences the absolute errors add, one must expect this procedure to be 
also ill-conditioned. If both Gq(z) and G(z) are known exactly, the difference Gq{z)~ 1 —G{z)~ 1 
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always gives a negative imaginary part for the self-energy as there would be a pole in G(z)^ 1 
for every pole in Go(z)^ 1 at the same energy with equal or larger residue. This is no longer 
guaranteed as soon as both Gq(z) and G(z) are only known approximately, and one has to use 
rather large values of the broadening parameter b to avoid unphysical oscillations in H u (z). 
This broadening in turn leads to a strong suppression of the high energy peaks because spectral 
weight is shifted from the center of the peak to its tails (to the high energy side due to eq. 

For the calculation of TF \z) via eq. (|]) on the other hand we do not expect to face these 
kind of problems that severely. Again, both quantities are calculated on the same basis by 



broadening the NRG results with (|10l) , i.e. with the same systematic error. This time, however, 
we divide them by each other, which means that only the relative errors will be propagated, 
leading to a numerically much more stable procedure. 

Let us support this rather qualitative argument in favour of expressing the self-energy as the 
ratio UF(z)/G(z) by comparing the spectral function A(u) obtained directly from the NRG 
(solid line in Fig. |Xp and the A(u) calculated from eq. (Q) with the self-energy expressed as in 
eq. (]8|) (dashed line in Fig. [1]). The spectral functions are defined as A(uj) = — ^QmG(uj + i0 + ) . 
The parameters are £f = — 0.1-D, U = 0.2D and -3mA(w + i0 + ) = Aq = 0.015-D, where 2D is 
the conduction electron bandwidth. For convenience we use D = 1 as energy scale for the rest 
of the paper. The differences between the two methods can be summarized as follows: 

• We find for the total spectral weight / <koA(uj) = w = 0.93 and / duA(u) = w = 0.9993. 
The 7% deviation in w can in principle be reduced by improving the resolution of the 
NRG calculation (smaller deviations have been achieved e.g. in ^ 0)- This is, however, 



not necessary in our case because the self-energy resulting from eq. (§) is an analytic 
function and the sum-rule w — 1 is then automatically fulfilled (apart from the very small 
numerical error). 

• The charge fluctuation peaks near E{ are much more pronounced in A(u). That the 
high energy features are usually underrated is a well-known problem in the calculation 
of dynamical properties with the NRG. This problem is at least partially resolved in our 
new scheme, since the main contribution in this part of the spectrum rather comes from 
the hybridization part A(z), which is treated exactly. 

• The oscillations in A(u) near uj = are due to the choice of the broadening which is 
obviously too small to see the correct low-frequency behaviour of the spectral function. 
These oscillations almost vanish in A(u) and the u 2 - behaviour can be clearly identified. 

• The deviation from the Friedel sum-rule 

M0) = - c * n+ , =:-jk (11) 

7r<577lA(«0 + ) 7tA 

(f» 21 for the parameters used here) is about 7% in A(u) and 4% in A(u). 

Although the error in the Friedel sum-rule is visibly reduced, the deviation is still a few percent. 
Its origin will be discussed in the following. 
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Figure 1: Impurity spectral function for £{ = — 0.1, £7 = 0.2, T = and — 3mA(w + i0 + ) = 0.015 
in units of the conduction electron band width. A(u) (solid line) is the result obtained directly 
from the NRG and A(u) (dashed line) is calculated via the self-energy eq. (§). The inset shows 
the region around the Fermi-level. 

Numerical aspects 

It is important to understand the origin of the deviation of A(0) (Friedel sum-rule) from its 
exact value as this lies at the heart of the numerical procedure. 

Typical results for A(u) and B(u) as calculated with the NRG are shown in figure |2|. Both 
spectral functions A(uj) and B(u) display a sharp resonance close to the Fermi energy. However, 
in contrast to A(u), which is positive definite and perfectly symmetric to u — due to the 
particle- hole symmetry, the function B(u) is obviously not positive definite and appears to be 
extremely asymmetric. 

As next step we must calculate the real parts, which are obtained via standard Kramers- 
Kronig transformation. The self-energy finally is given with eq. (§) as 

^ rj. ^ , 7 , TT ^eF(u + iO + ) +i%mF(uj + i0+) , . 

®,e£F(uj + i0 + ) + iQmE u (u + i0 + ) = U-—± -{ — — £ 12 

In the particle-hole symmetric case, $teG(ui + i0 + ) necessarily vanishes at u = 0, therefore 
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Figure 2: Spectral functions A(u) (solid line) and B(ui) (dotted line) for the same parameters 
as in Fig. 1 (directly from the NRG, not via the self-energy). The inset shows the region around 
the Fermi-level. 



which of course has to give the Hartree term U/2, and 



QmL u (i0 + ) = -U - ' ' (14) 

In the case of the standard single-impurity Anderson model we furthermore know that the 
Friedel sum-rule QmH u (i0 + ) = has to be fulfilled, which implies that 

oo 

KeF(i0+) = - / dujB(u)V- = , (15) 

J LU 

— oo 

where V(. . .) denotes the principal value. The relation (O) is obviously not trivial regarding 
the unusual shape of B(uj). It indeed turns out that 9ffeF(i0 + ) is numerically zero as long as 
the full spectrum of the Hamiltonian can be used. However, as soon as a truncation of states 
sets in the calculated value for 3?eF(i0 + ) suddenly jumps to a finite value, eventually leading 
to a violation of the Friedel sum-rule as observed e.g. in Figure [I]. 

This observation suggests that also high energy states are important to guarantee that 
?R,eF(i0 + ) = and that a slight violation of the Friedel sum-rule is almost unavoidable in this 
method. 
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3 Results 



Single impurity Anderson model 

As a simple example let us discuss the standard case of a constant 3mA(w + i0~ 



J A : 




I : 





(16) 



The application of the NRG to this model has been discussed very extensively in the literature 
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Figure 3: Real and imaginary part of the self-energy for E{ = —0.1, U = 0.2, and a constant 
A = 0.015. The inset shows the region around the Fermi level where the Hartree term was 
subtracted off the real part. 



||. [L0|. Thus the results presented here surely give no new insight into the physics of this model. 
They are mainly intended to give the reader a feeling for the quality of our method. 

Figure |^ shows the results for the real and imaginary part of ^ u (z) for £f = —0.1, [7 = 0.2, 
A = 0.015 and T = 0. As a first important point we note that the real part of the self-energy 
has a constant contribution. If we calculate $t,eYF(u) + 9fteS c/ (— lu), we obtain the expected 
value U/2 to within numerical precision for all u. This result also shows that, although B(u) is 
asymmetric, the final result obeys the particle-hole symmetry to a high precision. In addition 



the slope d^eTP \u) / 'du 
The imaginary part o 

as uj — > 0. In the vicintiy of the Fermi level we find the Fermi liquid property ^mYP (uo + i0 + ) oc 



uj=0 



is negative and large, corresponding to a high effective mass. 
(z) shows two pronounced peaks at u> ~ ±0.03 and a steep decrease 

Ui 
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uj 2 (inset of Figure |3]). However, as pointed out in the previous section, the Friedel sum rule 
QmS^ (i0 + ) = is not exactly fulfilled. The shift of SmS* 7 (i0 + ) ~ —0.0007 corresponds to a 
4% error in A(0) = 1/(ttA ). 

Figure [| finally shows the resulting spectral function for various values of U. As mentioned 




Figure 4: Spectral function for £f = —U/2, A = 0.015 and various values of U. The inset 
shows the resulting effective mass m* (filled circles) together with the expected behaviour 
m* oc VU exp(7rC//(8A )) (crosses). 

previously in section 2, we find pronounced charge fluctuation peaks at ±U/2 and the charac- 
teristic Abrikosov-Suhl resonance at the Fermi level. With increasing U this resonance becomes 
sharper. The corresponding energy scale expressed via the effective mass is shown in the inset 
to figure | together with the expected behaviour m* oc y/TI exp(irU / (8A )). 



Application to the Hubbard model 

The impurity Anderson model is not only useful to describe magnetic impurities in non- 
magnetic metals. It was shown only recently that in the limit of infinite spatial dimensions 
a lattice model (Hubbard model, Periodic Anderson Model, etc.) with local interactions can 
be mapped on an effective single impurity Anderson model. The quantity A (z), which in the 
single impurity model describes the coupling to the metallic host, becomes in general an en- 
ergy dependent quantity here, which has a meaning similar to the Weiss field in the mean-field 
theory of the Heisenberg model. Since A(z) is a dynamic quantity which must be determined 



self-consistently as functional of the one-particle self-energy [n], [T2|, [13]], the name "Dynamical 
Mean Field Theory" (DMFT) has been coined. 
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This self-consistency makes it necessary to calculate the self-energy E u (z) as accurate as 
possible. Here we want to demonstrate that the NRG together with the method of calculating 
S f/ (z) presented in the previous section is indeed a reliable and accurate method to do this job 
at T = 0. 

The first step in order to apply the NRG is the mapping of the impurity model on a semi- 
infinite chain for the case of a non-constant 3mA(w + i0 + ) which we have already described 
in fl5 |. As the resulting 3mA(w + i0 + ) can develop very narrow structures at the Fermi level, 
we need a reliable numerical method to calculate ~ 60 - 100 hopping matrix elements of the 
chain. This is done using arbitrary-precision fortran routines. Apart from the difference in the 
hopping matrix elements, the calculation of F(z), G(z) and ^ u {z) follows the same procedure 
as in the flat-band case. 



The simplest model for correlation effects in solids is the well-known Hubbard model [16 
This model is believed to have a rich phase diagram despite its comparatively simple form. 
DMFT studies at finite temperatures indeed revealed for example antiferromagnetic [|17| , [13 



and ferromagnetic transitions [jTj], W , Mott-Hubbard metal insulator transition [|T3|] etc. . Nev- 



ertheless, there still remain lots of interesting open questions, especially about the properties 
of the model at extremely low temperatures both at and away from half filling. 

Here we study the Hubbard model at T = for a semi-circular density of states po(e) 
corresponding to the Bethe lattice with infinite coordination number 

2 



p (e) = -^[Y^7\ (17) 



7T 



(JD — 1) at particle-hole symmetry and in the paramagnetic regime. The resulting spectral 
functions for the paramagnetic Hubbard model for various values of U are collected in Figure 
H With increasing U, the one-particle spectrum develops the typical three-peak structure with 
a quasiparticle peak at u = and the two Hubbard bands at ±U/2. Above a certain value 
U c ~ 2.93 the central peak vanishes and the system becomes insulating. 

Figure ^| and [F] show the real and imaginary part of the self-energy for the same parameters 
as in Fig. [5] (U — 1 and U = 4 are not shown). The Hartree term in the real part (= U/2) 
is subtracted. The negative slope at the Fermi level diverges as U — > U c . For U > 2.93 
(the insulating solution) the real part shows a l/cj-divergence. The corresponding 5-peak in 
the imaginary part is not plotted in Fig. [7| This 5-peak in lmH u (uj) emerges from a two- 
peak structure in the metallic regime, with the positions of the two peaks approaching u = 
for U — > U c . The imaginary part shows the Fermi liquid behaviour ImS f/ (u;) cx uj 2 at low 
frequencies for U <U C 

In order to give the reader an idea of the complex structures arising in lattice models, we 
show in figure |8| a comparison of the NRG flow diagram for the energy levels for the single 
impurity Anderson model with flat 3mA(w) (figure |8]a) and typical results for the Hubbard 
model in the paramagnetic metallic phase (figure |8|b) and paramagnetic insulating phase (figure 
^p). In contrast to the single impurity case the flow diagrams for the Hubbard model show 
a complicated crossover behaviour for high energies (low NRG iteration number) before they 
saturate into a fixed-point spectrum for large NRG iteration numbers, i.e. low energies. While 
these fixed-point spectra for the impurity model and the metallic solution of the Hubbard model 
(figures |8|a and b) are identical, i.e. both are Fermi liquids, the fixed point spectrum for the 
insulator (figure ||c) has a quite different structure (see eg. the flow of the first excited state 
with Q = 1, S = in figure [8]b and figure [8]c, Q is defined as the particle number with respect 
to half-filling) . In addition the behaviour of the hopping matrix elements for the three cases is 
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Figure 5: Local spectral function of the Hubbard model for various values of U. A quasiparticle 
peak develops for increasing values of U which vanishes at a critical value [7 C ~2.93, signalling 
the metal-insulator transition. 
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Figure 6: Real part of the self-energy for the Hubbard model (same parameters as in Fig. [5]). 
The negative slope at u = diverges at the metal-insulator transition. For U >U C , the real part 
shows a l/c<>divergence. 
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Figure 7: Imaginary part of the self-energy for the Hubbard model (same parameters as in Fig. 
5). A 5-function develops for U — > U c . 



shown in figure ||d (diamonds, circles and crosses, respectively). Note the oscillatory behaviour 
in the latter case characteristic for a system with a pseudogap. 



The data shown here are quite similar to those obtained by Georges et al. |T^] in that the 
quasiparticle peak seems to be isolated from the two Hubbard bands near U c . However, we 
always find finite spectral weight in the region between the quasiparticle peak and the Hubbard 
bands. For U >U C the quasiparticle peak vanishes but we do not see a real gap, i.e. a region 
with exactly vanishing density of states. The fact that there is no gap even above the metal- 
insulator transition is also visible in 5smH u (uj), as there is no region with vanishing 5smH u (u}). 
The strong suppression of the spectral density between the Hubbard bands and the quasiparticle 
peak for [7 = 2.85 mainly comes from the large values of dteH u (u). The question, whether a 
real gap will eventually emerge for higher values of U is currently investigated and a more 
detailed analysis of the metal-insulator transition at T = will be presented in a subsequent 
publication. 

For the time being we define the point where the transition from a metal to an insulator 
takes place by the divergence of the effective mass 



m* = 1 - ^-toe?F(uj) 

OUJ 



(18) 



of the quasi-particles. Note that this scenario completely neglects the possibility of a discon- 
tinuous transition for a U < U c . 

The behaviour of the effective mass as function of U is collected in figure |^. m* diverges at 
C/ c ~2.93 and the critical behaviour close to U c is consistent with a powerlaw with an exponent 
of ~ —2. Unfortunately, the data currently available do not allow a precise evaluation of this 
exponent. Note that our value for U c is considerably smaller than the value of U c % = 3.3 
mentioned in the work of Georges et al. ]13| . 
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Figure 8: Flow diagrams for the lowest energy levels Eg t g as function of the NRG-iterations N. 
The solid lines correspond to quantum numbers Q = 0, S =1/2 and the dashed lines to quantum 
numbers Q = l, S = 0. (a) The flat-band case with £f = — 0.2, U = 0.4, and a constant A = 0.015. 
For large N, the system flows to the Fermi liquid fixed point, while in the intermediate regime 
(N ~ 20) it is near the so-called local moment fixed point, (b) The Hubbard model with U = 2 
flows to the same Fermi liquid fixed point as in the flat band case, (c) The Hubbard model 
with {7 = 4 flows to the local moment fixed point corresponding to the insulating behaviour, 
(d) The hopping matrix elements of the semi- infinite NRG chain for Fig. ||a (diamonds), 
Fig. [8]b (circles) and Fig. |8|c (crosses), respectively. 

4 Summary 

In this paper we have presented a new method of calculating the self-energy of the single 
impurity Anderson model with the Numerical Renormalization Group method. In contrast to 
the standard approach where one calculates the self-energy from the Green's function alone, we 
express £ (z) as a ratio of two correlation functions. The central aspect of this paper is that 
this method is much more accurate than the usual method. 

The importance of this gain in accuracy goes beyond the mere improvement of the results 
for the single impurity Anderson model. Our method now in addition allows to apply the NRG 
to various lattice models within the Dynamical Mean Field Theory, where the self-energy of 
an effective impurity Anderson model has to be calculated self-consistently. As examples we 
recapitulated typical results for the single impurity Anderson model and presented results for 
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Figure 9: [/-dependence of the effective mass m* for the Hubbard model, m* diverges at 
U c ~2.93 which defines the critical value of the metal-insulator transition. 

the Hubbard model with a semi-circular density of states (corresponding to the Bethe lattice 
with infinite coordination number) at particle-hole symmetry and T = 0. 

As an interesting but not yet fully confirmed result in the latter case we find that the metal- 
insulator transition in this model appears more like a metal-semimetal transition, as there is 
no region in the spectrum where the spectral density exactly vanishes. We defined the critical 
U via the divergence of the effective mass and found a value U c ~ 2.93. 

A more detailed analysis of the metal-insulator transition, results for the Hubbard model 
away from particle-hole symmetry and the investigation of more complicated models (Periodic 
Anderson model, three-band Hubbard model, etc.) will be discussed in forthcoming publica- 
tions. 

We wish to thank T. Costi, J. Keller and D. Logan for a number of stimulating discussions. 
One of us (R.B.) was supported by a grant from the Deutsche Forschungsgemeinschaft, grant 
No. Bu965-1/1. We are also grateful to the EPSRC for the support of a research grant (No. 
GR/J85349). 

A The correlation function F(z) 

Here we present some details of the calculation of F(z) and discuss some of its properties. 

The NRG uses a discretized version of the Anderson model in a semi-infinite chain form (for 
details see J7|, §[). The resulting spectral functions will therefore be given as a set of discrete 
5-peaks. 
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The spectral representation of F(z) is 



5 H = \ E < n \flf\f\\ m >< m\f'l\n > 



x8 (u - (E m - E n )) ( 



g /3-En _|_ g — PE„ 



(19) 



The matrix elements < n\fifl f^\m >, < m|/||n > and the energies E n ,E m are calculated 
iteratively in the NRG. The two operators 



V 1 



1/2 
/2 



T/l/2 
K -l/2 



flfjfh 



(20) 



transform as 



1/2 



-?(g±l) 



(21) 



(q = ±1/2), with the spin operators 



fjfl, s- = f}f h 
: -x (/|/t - • 



(22) 



This allows us to use the Wigner-Eckart Theorem 



(Q, S, S z , w\ V q 1 ' 2 \Q', S', S'„ w') = (Q, S, w\\V^ 2 \\Q', S', w'){S', S' z , ± q\S, S z ). (23) 

The (Q, S, wlVq^WQ' , S', w') are reduced matrix elements and the (5", S' z , ~, q\S, S z ) Clebsch- 
Gordan coefficients. It is important to note that the operators V^ 2 transform in exactly the 
same way as the two operators 



VV l/2 

w 1/2 



P 



fl 



(24) 



This has the consequence that all the recursion formulas for the reduced matrix elements of 
Wg' 2 can be used for the calculation of the reduced matrix elements of V^ 2 . The only changes 
are in the particle numbers Q of the states involved and the initial values (see below). 

The states \n> and \m> in eq. (^) are classified in terms of charge Q (the total particle 
number relative to the half-filled case), total spin S, ^-component of the total spin S z and an 
additional label w 



\n > — \Qn, S n , S z>n , w n >, 



(25) 
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The sum over S Z)U and S Zym in eq. fll9|) can be performed exactly and we find 

B(u) = — J! (Q;^Wn||V r 1 1 / / 2 2 ||Q + l,£ m ,U> m ) 

Q,8,Wn 

x (Q + l, S m , wJflWQ, S, w n )5 (u - (E m - E g )) 
x ( e -^ +e -^)i=v / 25TI 

X | ^ ' S™, = S — 7} ^g) 

[ — V 5+1 : 5 m = S + 2 

We are only interested in the limit of zero temperature where we have 



^ , - _ „ .., 



x (Q g + 1, S m ,w m \\fl\\Qg,Sg,Wg)5(U-(E m -Eg)) 



for positive frequencies and 



x 5*3, u> fl |/JflQ fl - 1, S n , w n )5 (u - (Eg - E n )) 

v/2 V 9 \ /VTT : S n = S 5 + ± 

for negative frequencies. The ground-state is labelled by \g >= \Q g , S g , S Z) g,w g >, the ground- 
state energy is E g and the partition function Z reduces to the ground-state degeneracy. 
To set up the iterative diagonalization of the reduced matrix elements 

1 /2 

\Qni Sn, w n\\Vi/2 WQm, S m , w m ) we first of all need the initial values for the uncoupled impurity. 
The only non-zero matrix element is 

(0,^1^11, 0> = -l, (29) 

in contrast to the two initial values 



(0,^1/111 -i,o> 



2' 

(l,0f/|f0,~> = -V2. (30) 

1/2 ii 

Apart from the difference in the initial values and the fact that Q m = Qn + 1 for the < \V 1 j 2 \ > 
matrix elements, the recursion relations for both reduced matrix elements are identical and 
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given by 



N (Q,SM\V^\\Q',S',w') N = Y,Y, U Qs(™>rp)U Q/ s,(w',r>pi) 

rr' pp'=l 

■ N{Q,S,r; P \V^\\Q\S\r'-p') N , 



(31) 



with p,p' G {1,2,3,4}. The Uq^s are the unitary matrices which diagonalize the Hamiltonian 
matrix in the subspace with ccharge Q and spin S. The reduced matrix elements on the right 
hand side of eq. fl5I|) are given by 



N (Q,S,r; llVy^ 



N {Q,S,r;2 
N (Q, S,r; 2 
N (Q, S,r; 3 

N (Q, S,r; 3 
n(Q, S,r; 2 
N (Q,S,r;3 
N (Q,S,r;4: 



\Q + l,S±^,r';l) N = N . 



1/2 



1/2 



1/2 



1/2 



1/2 



1/2 



1/2 



1/2 



V, 



1/2 



1/2 



1/2 



1/2 



1/2 



1/2 



\Q + hS+- 
|Q + l,5-~ 
\Q + 1,S+- 



\Q + l,S-~ 
|Q + l,5-~ 
\Q+l,S+~ 
\Q + 1,S±~ 



r';2) N 
r';3) N 
r';2) N 



^(Q + l,S,r\\V 1 1 / / 2 2 \\Q + 2,S±±r') N _ 1 
2V 2S 2 +1 S N 1 {Q ' S ~ 5' rlV ^ lQ + h 5 ' " Vl 



1 

2' 
1 

2' 



- v-i (Q, 5 - ^, t\V$\Q + 1,5-1, r Vi 



- v-i (Q, S + ^, r\V$lQ + 1,5 + 1, r')^ 



2VS 2 + 5 

25+1 7v_i 
1 



(Q,S+± r\V$\Q + 1,5,/)^ 



-— — (Q, 5 - 1, r|V$lQ + 1, 5, r %_i 

ZD + IjV-l Z 



r';4)^ = w „ 1 (Q-l,5,r||^||Q,5±-,r%_ 1 



The spectral function B{ui) obeys the sum-rule 

/oo 1 
dwB(w) = -£ e -^» < n |/t/ t |„ >^< fj f] >, 
-oo Zj „ 



(32) 



(33) 



which can be easily derived by integrating eq. (O) over a;. In the particle-hole symmetric case 
this gives 

duB(u) = -, (34) 



(35) 



where we also find the following relation between B(u) and A(u) 

B(w)+B(-u) = A(u). 
This can be directly obtained from eq. (|6|). 
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